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ABSTRACT 



When a radar pulse impinges upon a target, the resultant 
scattering process can be solved as a linear time-invariant 
(LTI) system problem. The system has a transfer function with 
poles and zeros. Previous work has shown that the poles are 
independent of the exciting waveform and target's aspect, but 
they are dependent on the target's structure and geometry. 
This thesis evaluates the resonance estimation performance of 
two signal processing techniques: the Kumaresan-Tufts 
algorithm and the Cadzow-Solomon algorithm. Improvements are 
made to the Cadzow-Solomon algorithm. Both algorithms are 
programmed using MATLAB. Test data used to evaluate these 
algorithms includes synthetic and integral equation generated 
signals, with and without additive noise, in addition to new 
experimental scattering data from a thin wire, aluminum 
spheres and scale model aircraft. 
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I. INTRODUCTION 



When a radar pulse impinges upon a target, the resultant 
scattering process can be solved as a linear time-invariant 
(LTI) system problem. The system has a transfer function with 
poles and zeros. This description is provided by the 
Singularity Expansion Method (SEM) developed by Dr. Carl Baum 
[Ref. 1], of the Air Force Weapons Laboratory at Kirtland Air 
Force Base. Baum's SEM describes the induced current and 
scattered field using damped sinusoidal waveforms which 
resonate with complex natural frequencies unique to the 
object. These frequencies are determined by the object's 
composition and structural geometry. Natural frequencies are 
also the complex poles of the transfer function. These poles 
are independent of the incident electromagnetic excitation, 
including aspect and polarization, as initially postulated by 
Moffatt and Mains [Ref. 2]. Morgan [Ref. 3] has shown that a 
target scattering response (following illumination) can be 
represented as a weighted expansion of complex natural 
resonances whose poles are independent of the incident 
electromagnetic excitation. The problem of classifying a 
radar target can be solved, in principle, by determining the 
pole positions of the target's response. 
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Although this idea is not new, early attempts to 
demonstrate the practicability of such an identification 
system have produced poor results due to the presence of noise 
in the system. Two separate signal processing algorithms have 
been developed by Kumaresan-Tufts [Ref. 4] and Cadzow- 
Solomon [Ref. 5] to locate target poles with a high 
degree of accuracy, in the presence of noise. The Kumaresan- 
Tufts algorithm is intended for purely auto-regressive (AR) 
signals, while the Cadzow-Solomon algorithm is intended for 
auto-regressive moving-average (ARMA) signals. This thesis 
evaluates the viability of these two signal processing 
techniques by using new experimental electromagnetic 
scattering data. The use of this method improves 
identification of natural resonances in noisy signals with the 
correct selection of problem parameters. 

A. THE PROBLEM 

The approach employed to classify radar targets is 
comprised of two steps, based on natural resonances. The 
first step locates the position of the poles of the targets of 
interest. Numerical analysis of integral equation techniques 
will derive the poles for simple targets, e.g., a thin wire. 
For more complicated targets such as aircraft, the poles must 
be extracted from actual measurements of the target's response 
to incident electromagnetic excitation. The information 



2 



collected can then be used to form a database for comparison 
to the data obtained in actual field use. 

The second step for target classification is the 
comparison of the data obtained in field measurements with the 
data contained in the database. The target is classified 
based on the closest data match. One method of accomplishing 
the classification is to use the same signal processing 
technique that was used to form the database initially. Speed 
is an important factor in the classification process, as the 
time required to achieve classification must be less than the 
time for the target to become a threat. Another, more 
efficient way to perform the comparison is to employ the use 
of annihilation filters, as proposed by Morgan and Dunavin 
[Ref. 6]. This approach employs energy cancellation 
based upon the location of the target's poles. For example, 
a total target response (including early-time and late-time) 
may be subjected to a bank of annihilation filters. The 
filter that corresponds to the proper target will exhibit, in 
theory, a response coincident with the driven portion of the 
target response while annihilating the signal during the late- 
time portion of the response. 

A system used for classification of radar targets would 
require a bank of annihilation filters corresponding to each 
class of target. The system could then determine the class of 
target, based on the filter producing the minimal amount of 
output energy. 



3 



B . BACKGROUND 

Transient electromagnetic scattering can be described by 
showing an incident field in free space illuminating a 
perfectly conductive finite-sized object. Figure 1 depicts 
induced current on a scattering body. The magnetic field 
integral equation (MFIE) 
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Figure 1 Transient electromagnetic scattering 

describes induced current on the surface of a scattering body 
as: 

J{r,t)=2fixH 1 (x,t)+f[ K(r,'r / I t) • J(i, t- l r ~ r I ) -els (D 

JJ Spy c 

where , 

• J is the electric current density, 

• i? is the unit normal vector outward to the surface. 
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• H 1 is the incident magnetic field intensity at the 
surface, 

• K is the dyadic Green's function kernel 

and, the principal-value (PV) type integral excludes the point 
"r^r 1 . The cross product, 2 fixH 1 , represents the physical 
optics portion of the induced current. The principal-value 
surface integral, S pv , represents the contribution due to 

induced physical optics currents other than the contribution 
from the cross product. This current represents "feedback" 
effects from currents at other points on the object. 

When H 1 = 0, Equation (1) solutions become the source-free 

("natural") current modes of the scattering problem. These 
modes may be described by the sum of the exponentially damped 
sinusoids, or: J^(r) exp (s n t) , where s n are the poles or 

natural resonance frequencies in complex conjugate pairs of 
the form 

< 2 > 



where , 

• o D = the damping rate in Nepers/sec, 

• o> n = the frequency in radians/sec. 

These resonance frequencies are functions of the geometry and 
composition of the scattering object. Although the number of 
poles is, in theory, infinite for any given object, only a 
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finite subset will be excited, due to the finite bandwidth of 
the incident field. 

Figure 2 represents a plane wave impulse illumination, in 
which part of the object has been illuminated from the 
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Figure 2 Plane wave impulse illumination 



incident field wavefront. The scattered field is composed of 
two parts: the early-time driven response and the late-time 
natural mode response. The early-time driven response occurs 
as the excitation wavefront travels over the object and each 
point becomes excited. The late-time natural mode response 
occurs, by definition, when the excitation has been completed 
over the complete body. Throughout the early-time phase, the 
impulsive plane wave incident field will be identically zero 
at all points on the surface except on the conformal ring 
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where the intersection with the wavefront occurs. This area 
is indicated by the dotted line in Figure 2. This conformal 
ring changes cross-sectional shape and position as the 
wavefront moves. The induced surface current on the ring is 
therefore composed of the physical optics term plus the 
feedback current, as described in Equation (1) . Thus, no 
induced current is present at the area of the object ahead of 
the wavefront. 

The back-scattered far-field, resulting from the surface 
current on the object, will be of the form 

s ' ( - rj5 ' c> shs-h /j>* 7(? - t - SL p- ) o) 

S 

The unit vector indicates the direction of the plane wave's 
propagation. The back-scattered far-field for a fixed point 
is then a result of integrating the current at every point on 
the object's surface. Thus, the back-scattered far-field will 
be of the form 

cm 

H s (-rj5, t) =u(t- — ) t) + £ H n {-rjS, t) exp (s n t)] . ( 4 ) 

C n— « 

n* 0 

Two individual terms emerge from these calculations. The 
first term in Equation (4) represents the physical optics 

scattered field generated by the 2 fixH 1 driven current. The 
second term represents the scattered field produced by the 
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source-free wake current or "feedback” current. The physical 
optics scattered field is highly aspect dependent. 

During the early-time portion of the target’s response, 
the scattered field produced by the "feedback" current term in 
Equation (4) contains exponential resonance terms with time- 

varying coefficients H n , as described by the Singularity 

Expansion Method (SEM) [Ref. 1]. This form of the SEM 
expansion is termed "Class 2" and is presented analytically by 
Morgan [Ref. 7]. For the monostatic scattering case, 
the transition from early-time to late-time in the target's 
response will occur at 

A t=x+ (1+cosa) •— , ^ 

c 

after the leading edge of the scattered pulse returns to the 
radar antenna, where, 

• t is the incident pulse width, 

• D is the target's largest dimension, and 

• a is the aspect of the target. 

The term a also represents the angle between the target's 
largest dimension and the direction of wave propagation, where 
c is the velocity of propagation or speed of light. At this 
transition time, the physical optics field vanishes. During 
the late-time portion of the target's response only the second 
term in Equation (4) remains. However, it is now comprised of 
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constant coefficients, H a . This form of the late-time SEM 

expansion is termed "Class 1". 

The early-time scattered field is therefore constituted of 
both a physical optics term and a "Class 2" SEM expansion with 
time-varying coefficients, while the late-time scattered field 
is described by a simple "Class 1" SEM expansion with constant 
coefficients. 
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II. POLE EXTRACTION ALGORITHMS 



Obtaining the target response is the first step to be 
completed in a Non-Cooperative Target Recognition (NCTR) 
system. Once the response is obtained in time domain, the 
target's poles may be located. Numerically, these poles may 
be described by the damping rate o and the radian frequency 
to or, in other words, the complex number indicating the pole's 
location in the s-plane or z-plane. The location of the poles 
may also be described analytically by solving the boundary- 
value electromagnetic problem. This may be easily done only 
for simple shape targets like thin wires and spheres. 
However, for general targets, detailed knowledge of the 
target's composition (dimensions and materials) and, also, 
access to a large amount of computing power is required. 

The next section of this thesis discusses the analytical 
method which is used to compare the poles extracted from 
measurements taken with a thin wire. The viability of the 
Cadzow-Solomon algorithm will also be tested. 



A. EARLY METHODS 

The NCTR system is used to discriminate and to classify 
targets with relatively similar shapes and dimensions. In 
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such cases, the location of the poles from these targets will 
be relatively close, requiring a high degree of precision in 
measurement and estimation. A basic method in signal 
processing is to obtain the signal's spectrum using the Fast 
Fourier Transform (FFT) . The FFT algorithm yields results in 
a short amount of time. However, the FFT may be used only as 
a tool to add information to the basic methodology when 
locating radar target natural resonances. This is because the 
frequency resolution, equal to the reciprocal of the time 
interval between the samples, is of the order of several MHz 
or higher. Another limitation of the FFT is that it provides 
only real frequencies within a signal, while both the damping 
rate and the frequency of the poles are required for target 
classification. 

Thus, other methods providing more accurate results needed 
to be developed. 

1. Direct Minimization 

As previously described, the transient scattered 
signal waveform will have the form 

y(t) =y £ ( t) [u( t) -u(t-T 0 ) ] +y L ( t) -u(t-T 0 ) +N( t) (6) 



where , 

• y B (t) is the early-time portion, 

• y L (t) is the late-time portion and 

• N(t) is the undesirable noise and clutter. 
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Morgan [Ref. 7] presents the way to determine the natural 
resonance parameters by modeling the late-time scattering 
response as a sum of damped sinusoids, 

«• 

=52 A n * e xp(o„t) -cos (o> n t+4> n ) (7) 

12-1 



where j? represents the modeled signal and A a ‘exp(j$ n ) 

represents the aspect dependent residues. The digital domain 
form of the model is 

N 

?(k’A t) 'Vexp (o D 'kA t) *cos t+<t> n ) (®) 

12-1 



where , 

• A n is the amplitude, 

• a D is the damping rate, 

• < a n is the radial frequency, and 

• <t> n is the phase. 

After comparing the modeled signal f n to the actual received 
discrete signal y B , the mean square error may be obtained as 

e n={y n -? D ) 2 . The four parameters of the model (A n , 4> n/ a n , u> n ) 

may then be adjusted in order to derive the mean square error 
minimum. As this is a multi-dimensional, highly non-linear 
minimization problem, it is computationally inefficient. 
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2. Prony's Method 

Blaricum and Mittra [Ref. 8] present a novel 
approach for systematically deriving the complex poles and 
residues of a target from a set of time domain data. The 
method is based on Prony's algorithm, used to model the late- 
time response of a radar target. The set of time domain data 
is the discrete set of sampled transient values of the impulse 
response J ( t n ) or I D . The method is based on the fact that 

I B must satisfy a difference equation of order N, written as 

E«,V4n' (5) 

p* 0 

where p+k=0 , 1, . . . , 2N-1 and 2N are the data samples used. This 
leads to the solution of a matrix equation 



J 1 


I 2 ... 














J 3 ... 


^1 
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a N - 1 
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*n* 2 




I N*1 — 






. 3 1 . 







After the coefficients a L are found, the method solves the 
roots of a polynomial as 

z N -a 1 z N ~ 1 a w =0 . (ll) 

If d m (m=l, 2 , . . . ,N) represents those roots, then the poles of 
the system model in the z-plane are the roots d m . The poles 
in the s-plane may be found using the formula 
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ln(c?J 

At ' 



( 12 ) 



s 



where At is the time-stepping interval used in obtaining the 
sampled data. The matrix in Equation (10) is in the form of 
a transposed Vandermonde matrix, whose inverse can be computed 
in closed form. This method requires at least 2N equally 
spaced transient data samples to find N poles. If greater 
than 2N samples are desired, then one may obtain a least 
squares type fit to the matrix in Equation (10) . Since the 
system order is not known a priori in the NCTR problem, 
Blaricum and Mittra [Ref. 9] present two schemes for 
determining the number of poles contained in the transient 
response. The first is the Householder orthogonal ization 
method, and the second is an eigenvalue method. When the data 
are noisy, Blaricum and Mittra [Ref. 9] indicate that the 
eigenvalue method is the better method. This method will be 
discussed in Section 3 of this chapter when a bias 
compensation method is examined. 

3. Kumaresan-Tufts Algorithm 

Kumaresan and Tufts [Ref. 4] modified Prony’s method 
by using the "backward linear prediction" technique and 
"singular value decomposition" (SVD) to alleviate the 
sensitivity to noise and the need for a priori knowledge of 
the system order. 
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Analytically, Kumaresan-Tufts algorithm modifies 
Prony's method as follows: 

• The system order is intentionally overestimated. The 
excess order provides the flexibility to the system to 
model the noise, improving the estimation of parameters of 
exponentially damped sinusoidal signals in noise. 

• Singular value decomposition alleviates severe ill- 
conditioning of the data matrix. 

• The causality of the system is used to separate the 
computed signal's poles from the spurious computed noise 
poles introduced by the overestimated system order. 

The separation of the signal poles from the noise 
poles is the result of the "backward prediction" that causes 
the signal poles to fall outside the unit circle in the z- 
plane, while the extraneous "noise poles" fall inside the unit 
circle. 

Moderately large values of system order are essential 
in improving the accuracy of the pole location estimates. 

a. System Model 

Suppose that the N samples of the observed time 
domain data of a response signal y(n) consists of samples of 
M exponentially damped signals in white Gaussian noise w(n) , 
as represented by 

M 

y (n) a**exp (Sj^n) +w(n) , n=0, l, . . . ,N-1 . (13) 

Jc-l 



The following linear prediction equations can be formed, using 
real time domain data in the backward direction. 
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( 14 ) 



' y (2) y(3) ... y(L-t-l)- 




a ( 1 ) ’ 




' y(l) ' 


y (3) y (4) ... y(L+2) 




a (2) 


= _ 


y( 2) 


y(N-L+ 1) y(N) . 




a(L). 




y(N-L)' 



The terms of the above equation may be also represented by 

D-a=-h (14a) 

where , 

• D is the data matrix (N-L) xL, 

• a is the vector of backward prediction coefficients 
(Lxl) , and 

• h the data vector (N-L)xi. 

In the above matrix equation, L represents the overestimated 
system order, chosen to satisfy the inequality M<L<N-M. The 
matrix equation is solved with the SVD method, using the 
pseudoinverse matrix D * , as a=D*-(-h) . Here the coefficients 
a ± correspond to those in Equation (11), where the roots of 

the polynomial define the poles of the system model in the z- 
plane. 

As with Prony's method, the Kumaresan-Tufts method 
is intended for purely auto-regressive (AR) signals. Both 
methods, therefore, use the late-time portion of the response 
signal when addressing the target classification problem. 
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b. Singular Value Decomposition 



Singular value decomposition (SVD) is the basic 
technique on which the Kumaresan-Tufts algorithm is based. 
There are two main advantages with SVD: 

• It helps to alleviate the effects of ill-conditioning of 
the data matrix. 

• It separates the signal poles from the noise poles. 

The following discussion of SVD is a synopsis of material 
taken from Strang [Ref. 10] . 

The SVD is closely associated with the eigenvalue- 
eigenvector factorization of a symmetric matrix : A=Q-A-Q t . 
SVD makes a similar factorization, but for any (mxn) matrix 
A, as follows: 

a=q 1 -T,-q 2 t < 15 > 

where Q lt Q 2 are two orthogonal matrices and S is the diagonal 

(but rectangular) matrix with its positive entries (also 
called sigma, in the form of o x , o 2 , . . . , o r ) . These entries are 

the singular values of A, filling the first r places on the 
main diagonal of E, where r is the rank of A. These entries 
sigma are also the square roots of the nonzero eigenvalues of 
both AA T and A T A . The columns of Q x (mxm) are eigenvectors of 

AA T , and the columns of Q 2 (nxn) are eigenvectors of A T A . 
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The SVD works well for numerically stable 
computations. An initial reason is that Q 1 and Q 2 are 

orthogonal matrices that never change the length of a vector. 
Since \\Qx\\ z =x T Q T Qx=\\x\\ 2 , multiplication by Q cannot destroy 
the scaling. A second reason is that SVD gives a more stable 
measure of the rank of A. 

The prime use of SVD is to solve every linear 
system in the form of Ax=b. For every matrix A (m x n) , which 

can be factored into A=Q 1 -H'Q * , another matrix can be defined. 
This matrix will be the pseudoinverse of A, as follows: 

A*=Q 2 , 'L*'Qi ( 16 ) 



where, Q x , Q 2 are the same orthogonal matrices found with the 



SVD. The value E* is the diagonal (n x m) matrix, with its 

The 



11 1 

entries on the main diagonal being — , — 



°1 °2 



optimal solution of Ax=b is 

x+=A+'b ( 17 ) 

that is, the minimum length solution to the nearest equation, 
Ax=p. As the column space of A* is in the row space of A, 
then x* is always in the row space of A. 
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To solve the system of equations in Equation (14) , 
the pseudoinverse of matrix D may be found as in Equation (16) 
in the form D*=V'H*’U T , where E + will be a (Lx (N-L) ) matrix, 
whose L singular values are the reciprocals of those found in 
the E matrix. The minimum length least squares solution of 
D'a=-h is a=D*'(-h) =V-'L*-U T -(-h) . 
c. Bias compensation 

In the discussion of the SVD above, S is a (N-L) xl 
diagonal matrix, with its entries being the square roots of 
the nonzero eigenvalues of both DD T and D T D. In the noiseless 
case, the prediction equations are satisfied exactly. If the 
system order has been overestimated as L, when M is the actual 
order of the system, the diagonal of E splits into a signal 
subspace with M positive singular values of o 1 ,o 2 , . . . ,o M and 

a subspace with L-M zero values. For the noisy signal case, 
the first M positive signal values are perturbed into a noisy 
signal subspace 

o 1 ,o 2 , . . . ,o M with eigenvalues X^X^. • . ^Xj* , ( 18 > 

and a noise subspace with L-M values 

°m*i • °m*z • • • • >°l eigenvalues X« +1 *X^ 2 * • • • 

Kumaresan and Tufts noticed [Ref. 4] that the noise 
perturbed the signal's singular values, the reason for the 
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bias towards the unit circle in the z-plane in the pole 
position estimates. Kumaresan and Tufts proposed a 
compensation method which moves the poles back towards the 
center of the z-plane. This method is based on averaging the 
smallest L-M singular values, o w+1 , °m* 2 > • • • • °l> an< * then 

subtracting this average value from each of the first M 
singular values, The smallest L-M values are 

then set to zero and a new matrix, 2 , is used to compute the 
pseudoinverse D * . Although Kumaresan and Tufts claim that 
this method gives better results, the testing completed for 
this thesis shows that the bias towards the center of the z- 
plane is greater than the bias error that the noise makes. 
This result causes much more perturbation on the pole position 
estimates. 

Norton [Ref. 11] proposed another 
compensation method based on the Eigenvalue Shifting Theorem. 
The noisy data matrix may be described by D=S+N, where S is 
the noiseless signal data matrix and N is the noise data 
matrix of the wide-sense stationary white noise process, v it 

written as 



v i 






N= : 



V 



N-L 



”N 



( 20 ) 
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The expected value of DD T may be found from 



DD t =E[ ( S+N) ( S+N) n =E[SS T ] +E[S’^ r ] +^[^5^ +E[NN T ] . (21) 

As white noise has a mean of zero, the two terms EiSN 7 ] and 
2?[-M£> r ] become zero. Since the noise is wide-sense stationary, 
ElNN 7 ^ =ol'I , where o v is the noise variance and I is the 

identity matrix. As S is deterministic, E[SS T ]=SS T . This 
leads to the formula 

E[DD T ]=SS T *al'I. (22) 

The eigenvalue shifting theorem [Ref. 11] describes the case 
when the eigenvalues of SS T are X d , where the eigenvalues of 

E[DD H are X d +al . Therefore, the eigenvalues of DD T , which 

are the squares of the singular values of D, are increased by 
the noise variance, oj . 

In the noiseless case, these singular values would 
be zero. If the noise singular values are squared and 

averaged, an unbiased estimation of oj may be obtained. 

Norton [Ref. 11] proposed correcting the signal singular 
values by first subtracting the noise variance oj from the 

squares of the first M diagonal entries of matrix S and then 
taking the square root of the reduced values as the new 
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diagonal entries of matrix S (while the noise singular values 
are set to zero) . As in the Kumaresan -Tufts bias compensation 
method, the pseudoinverse D* can be found from the new matrix 
E. 

Although Norton's method seems to be a more correct 
bias compensation method, both methods are based on the fact 
that a decision has to be made about the actual order of the 
system. The separation between the signal eigenvalues and the 
noise eigenvalues is not readily obvious, because the 
eigenvalues of matrix DD T or D T D appear to steadily decrease. 
d. Performance and Earlier Results 

Kumaresan and Tufts tested the algorithm, obtaining 
good results using synthetic data with various levels of white 
Gaussian noise, as low as 15 dB SNR. Norton [Ref. 11] 
developed the algorithm as a computer subroutine and tested it 
with various sets and types of data. Synthetic data were 
tested at various SNR values, ranging from 90 dB to 7 dB. 
When the SNR decreased, the error in pole position estimates 
increased. Norton claimed that the algorithm gave good 
results for SNR>20 dB. Norton also used simulations of thin 
wires produced by Morgan's Time-Domain thin wire integral 
equation (TDIE) computer program to test the algorithm. The 
results of those tests illustrated the aspect independence on 
estimated poles of the target, [Ref. 11]. 
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Larison [Ref. 12] programmed the algorithm 
using Fortran and tested synthetic and thin wire integral 
equation data at various SNR's ranging from 90 dB to 7 dB. 
Like Norton, he claimed good results in extracting the low 
frequency poles position with SNR as low as 20 dB by using 
Norton's bias compensation method. Larison maintained that 
the two most critical parameters are: 

• to select the appropriate starting point to begin the data 
processing 

• to select the appropriate system order to provide the best 
possible results. 

The algorithm has been tested by the author of this 
thesis, after developing the algorithm using MATLAB. The 
algorithm was tested using synthetic data at various SNR's 
ranging from 30 dB to 10 dB, without using the bias 
compensation method. The synthetic signal response is based 
on ten pole pairs within a frequency range of 1-10 GHz, with 
a medium Q damping factor using k=0.7. This chosen value of 
k simulates typical expected levels of damping from measured 
scattering responses of real targets (e.g., the thin wire). 
Damping factor and SNR level are discussed, in detail, in 
Chapter III, Section A, of this thesis. Table I shows the 
poles and residues used for the testing. 

Figures 3 to 10 illustrate the evaluation of the 
Kumaresan-Tufts algorithm. Figures 3 to 6 illustrate the 
synthetic signal waveform for various SNR's. Figures 7 and 8 
illustrate the spectrum of the synthetic signal for SNR=30 dB 



23 



and SNR=10 dB. Figures 9 and 10 illustrate the true and the 
algorithm poles in both the s-plane and z-plane, for various 
SNR's. For the above simulation, the values used for N and L 
were 300 and 60, respectively, while M (actual order) was 20. 
The simulation revealed that as the noise increases relative 
to the signal, there is a bias of all the poles towards the 
unit circle, specially the higher frequency poles. 



Table I MEDIUM Q SYNTHETIC POLES 



f n 


a n 


w n 


A n 


*» 


[GHz] 


[GNp/s] 


[Grad/s] 






1 


-0.3562 


6.2752 


1 


0 


7 


-0.7124 


12.5504 


1 


0 


S 


-1.0687 


18.8256 


1 


0 


1 


-1.4249 


25.1007 


1 


0 


5 


-1.7811 


31.3759 


1 


0 


6 


-2.1373 


37.6511 


1 


0 


7 


-2.4935 


43.9263 


1 


0 


8 


-2.8498 


50.2015 


1 


0 


9 


-3.2060 


56.4767 


1 


0 


10 


-3.5622 


62.7518 


1 


0 
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Figure 4 Generated synthetic signal with 30 dB SNR 
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Figure 5 Generated synthetic signal with 20 dB SNR 
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B. CADZOW- SOLOMON ALGORITHM 

Many investigators have used both excitation and response 
data to identify linear systems. Cadzow and Solomon [Ref* 5] 
proposed an identification procedure in which both the 
excitation and response are contaminated by white noise. This 
method is based upon the null space characterization of an 
associated "data matrix". The Cadzow-Solomon algorithm is 
used in this thesis to simulate and demonstrate the 
identification problem. 

1. System Model 

The Cadzow-Solomon algorithm is based upon the Auto- 
Regressive Moving-Average (ARMA) model. In an ideal modeling 
situation, the assumption may be made that the excitation 
signal x(n) and the response signal y(n) are perfectly 
related by means of an ARMA relationship, such that 

L K 

yin) =£ a^in-i) +Vb i x(n-i) , (23) 

i- l i^o 



as associated with the transfer function 



H(z) 



b 0 +b 1 z~ 1 +. 



. +b K z~ K 



(24) 



where , 

• a i are the coefficients which correspond to the poles of 
the transfer function, and 
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• b ± are the coefficients which correspond to the zeros of 
the transfer function. 

The order of the numerator of the system's transfer function 
is K and the order of the denominator is L. The classic 
modeling problem is to identify the system's coefficients, a. i 

and b it from a finite set of observations of the excitation, 

x{n ) , and response, y(n) , time series. 

Cadzow-Solomon present the set of equations for the 
algorithm in matrix form as follows: 



'y(o) y (-1) ... y(-p) ' 
y(l) y(0) ... y(l-p) 




1' 




'x(O) x ( —1 ) ... x(-g) 

x(l) x(0) ... x(l-g) 






y(N) y(N-l) ... y(N-p ) , 




a p. 




x(N) x(N-l) ... x(N-g) 




w 



or, more concisely 



where , 

• X q is the (N+l)x(q+i) "excitation data matrix", 

• Yp is the (N+l)x(p+l) "response data matrix", 

• a p is the (p+l)xi "auto-regressive parameter vector" and 

• b q is the (q+l)xi "moving-average parameter vector". 

The least squares solution for a causal, real data, 
overestimated ARMA model (p*p,g*g) can be found by taking a 
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linear algebraic approach based on the eigenvalue-eigenvector 
decomposition of the data matrix D p q = [ Y p i -X Q ] as 

Dp, qD Pt q 'U k -X k 'U k , lzkzp+q+2 . ( 27 ) 



The parameter vector is found from the procedure to be 




1^(1) l 2 



u k ( 1) _ 

E -nr - * 17 * 

*•1 A k 



( 28 ) 



where u k is the eigenvector associated with the eigenvalue X*. 

If the actual system order is (p, g) , then the system order is 
less than the overestimated order (p, q) . At least 
s=l+min (p-p, q-q) of the eigenvalues in the decomposition, as 
described in Equation (27) , must be identically zero for the 
noise-free case. 

Following the technique of "backward prediction" as 
used in the Kumaresan-Tufts algorithm, the matrix equation 
(Equation 25) may be modified as follows: 

y (2) ... y (L+l) 

y (3) ... y (L+2) 

y(N-L+ 1) ... y(N) 



x(l) 

x(2) 



x(K+l) 
x(K+ 2) 



x(N-L) ... x(N-L+K ) 






y(D 

y(2) 

y(N-L) 



( 29 ) 



or, in matrix notation 







( 30 ) 
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This matrix is associated with the coefficients and b ± as 

shown in Equations (23) and (24) . The size of the excitation 
and response data matrices are (N-L)x(K+l) and (N-L) xl, 
respectively, while the size of the ARMA parameter vector is 
(L+K+l) xl . 

As in the Kumaresan-Tufts case, singular value 
decomposition can also be used in Cadzow-Solomon's algorithm 
to provide the minimum-norm solution. Its use reduces ill- 
conditioning of the data matrix, D yx , and separates the signal 

poles from the noise poles across the unit circle. The 
optimal solution of Equation (30) is 



a 



b\ 



[ D yxY'y- 



(31) 



2 . Bias Compensation 

When the system's order, (p, g) , exceeds the true order 
of a noise contaminated system, (p, g ) , for p>p and g>g, 
Cadzow and Solomon maintain that there will be 

s=l+min (p-p, g-g) ( 32 ) 

eigenvalues, which will asymptotically approach (N+l) • o* for 
large N, where o* is the noise variance. The eigenvalues in 
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the noise-free case are zero, and therefore, the parameter 
vector may be found from the equation 




a 



(33) 



where the u k terms are the eigenvectors associated with the 
smallest (s) eigenvalues, X k . 



as described above because 

• the excitation noise, o v , is zero for the radar target 
classification problem, and 

• the g-g, cannot be directly determined since the 
extraneous zeros cannot be separated from the signal zeros 
in the same manner as the poles are. 



algorithm is that additive noise is different for input and 
output data. Norton [Ref. 11] described that if the input 
data noise is v i and the output data noise is w A , the noisy 

data matrix may be expressed as 



However, the above method may not be applied exactly 



Another problem occuring in the Cadzow-Solomon's 



Dyx \Py '• -Px] Sy X + N yx 



( 34 ) 



where, N yx =[N y ; N x ] , or in matrix form 



*y= 



V, 



W N-L * 1 



W, 



L* X 



V, 



N 



( 36 ) 



Therefore, the correlation product will be 



Md.... dZ\ =s._ sZ+Mn.. 



( 37 ) 



As the additive noise is different for the input and output 
data, the noise correlation matrix ElN^N^] is not of the form 

o 2 T. So, the eigenvalue compensation method that Norton 

described as being applicable in the Kumaresan-Tufts algorithm 
is not applicable in the Cadzow-Solomon algorithm. 

By examining the diagonal entries of the data matrix 

product, Dp tq 'D pq , it may be seen that the first (g) entries 

are close to 1. In the noise-free case, these first (g) 
diagonal entries are exactly equal to 1. This is a direct 

consequence of the fact that the diagonal entries of Dp iq 'D pq 
contain bias errors proportional to the noise variances a 2 and 

al. Since the excitation data used in radar target 

classification problems are noise-free, these bias errors are 
proportional only to a 2 w . The bias compensation may be 

obtained by setting the first (g) diagonal entries equal to 1 
and then finding the eigenvalues and eigenvectors from the 



34 



The parameter vector may then be 



corrected matrix, Dp >q 'D p q . 
found from 



a 



p 




p f . 2 IMP I* 

Jc-l ^-Jlr 



-1 



P*<J + 2 

E 

k*l 




(38) 



3. Earlier Results 

Norton [Ref. 11] programmed the Cadzow-Solomon 
algorithm in Pascal and tested it on various types of data. 
Using synthetic data, Norton tested the algorithm using three 
SNR values, 30 dB, 20 dB, and 10 dB. The results were good at 
the 30 dB level, but as the SNR decreased, the error in the 
pole position estimates increased. Norton also used 
simulations of thin wire scattering produced by Morgan's TDIE 
computer program to test the algorithm, obtaining better 
results than those obtained from the Kumaresan-Tufts 
algorithm. These results occured because the Cadzow-Solomon 
algorithm takes into account only the input and output data 
and does not require purely late-time data. 

Larison [Ref. 12] programmed the Cadzow-Solomon 
algorithm in Fortran and tested it using synthetic data and 
TDIE data, at various SNR's, ranging from 90 dB to 7 dB. 
Using bias compensation, Larison claimed good results in 
extracting the low frequency pole positions with SNR as low as 
20 dB. 
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Murphy [Ref. 13] tested the Cadzow-Solomon 
algorithm using Larison's Fortran programs, using synthetic 
data as well as thin wire integral equation generated data and 
measured data. Murphy generated three separate signals, each 
based on ten pole pairs covering the frequency range of 1-10 
GHz. Depending on the level of damping for each signal, 
Murphy obtained an average error between the true and the 
extracted poles having values on the order of 10* 2 for the 
values 

• SNR=10 dB and for HIGH Q, 

• SNR=20 dB and for MEDIUM Q, and 

• SNR=30 dB and for LOW Q. 

Murphy made the observations that, 

• The best results were obtained by choosing a starting 
point located within several points of the zero crossing 
nearest to the first obvious response of the excitation. 

• The best results were obtained by using a data matrix in 
which the overestimated number of poles to the true number 
of poles was of the order of 2.5, while the number of 
asking zeros was equal to the number of true zeros. 

• The best results were not always obtained when using the 
bias compensation scheme as proposed by Norton. 

When processing the thin wire data, Murphy calculated the 

feed-forward order of the system by determining the length of 

early-time as equal to 2L/c. Recent research at NPS, that is 

as yet unpublished, indicates that this method is' not correct, 

as each excited point along the target excites its adjacent 

point. This early-time may be described as (1+cosa) • L/c, as 

shown in Chapter III, where a is the aspect of the target. 
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When processing the TDIE data, Murphy did not obtain good 
results for SNR=20 dB, in spite of the fact that the poles 
extracted from thin wire measured responses appeared to give 
good results for the low frequency poles. 
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III. ALGORITHM TESTING 



The initial objective of this thesis was to evaluate the 
viability of the Kumaresan-Tufts and the Cadzow-Soloinon 
algorithms. In a radar target classification problem the 
finite duration excitation signal produces both early-time and 
late-time response signals. In this situation, the Cadzow- 
Solomon algorithm is more significant. Thus, the main effort 
of this thesis was changed to evaluate and improve the Cadzow- 
Solomon algorithm using both synthetic and real data. The 
Kumaresan-Tufts algorithm was evaluated and tested only with 
synthetic data, as presented in Chapter II. 

The Cadzow-Solomon algorithm was tested in two phases. 
The first phase of testing used synthetic data, while the 
second phase was performed with thin wire measurement data. 
The synthetic data testing phase attempted to simulate the 
conditions expected from the response of a simple target 
during the presence of a stationary white noise. The thin 
wire data testing phase attempted to evaluate the conditions 
appearing from a real target response. Those conditions were 
then compared with those simulated from the computed Time 
Domain Integral Equation (TDIE) thin wire response. 
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A. SYNTHETIC SIGNAL MODEL 



As the Cadzow-Solomon algorithm is based on the ARMA 
model, the representation in Equation (23) has been used to 
produce the synthetic signal response. 

This signal response is based on ten pole pairs within a 
frequency range of 1-10 GHz, with a medium Q damping factor 
using k=0.7. The poles were developed in accordance with the 
following equation: 



-^--^-ln(Jc) 

2n 



( 39 ) 



Table II lists the s-plane poles used in the synthetic signal. 



Table II MEDIUM Q SYNTHETIC POLES 



f n 


° n 


w n 


[GHZ] 


[GNp/s] 


[Grad/s] 


1 


-0.3562 


6.2752 


2 


-0.7124 


12.5504 


3 


-1.0687 


18.8256 


4 


-1.4249 


25.1007 


5 


-1.7811 


31.3759 


6 


-2.1373 


37.6511 


7 


-2.4935 


43.9263 


8 


-2.8498 


50.2015 


9 


-3.2060 


56.4767 


10 


-3.5622 


62.7518 



The chosen value of k simulates typical expected levels of 
damping from measured scattering responses of real targets, 
(e.g., the thin wire and scale model aircraft targets). The 
sampling frequency used to convert the s-plane poles to z- 
plane poles was 51 GHz, based on N=256 samples over a time 
window of t 0 =5 nsec: 
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1 N-l 



(40) 



fs= 



At 



to 



1. Coefficient Generator and Recursive Signal Generator 
The coefficients a i are obtained by multiplying the 

terms (z-z x ) ( z-z 2 ) . . . ( z-z 20 ) , where z i is the i th z-plane pole 
associated with the s-plane pole in the relationship 

z i =exp [ (Oj+jQj) -A t] , (41) 



and equal to the product of those terms with the polynomial 



z 20 -a 1 z 19 -a 2 z 18 - . . . -a 20 , 



(42) 



The coefficients b i are obtained by the inverse partial 



in- 



fraction expansion using the term — , where Rj, is the i 



th 



Z- Zi 



residue of amplitude value 1 and phase difference 0 for each 
of the signal poles. The time domain signal response of the 
ARMA model is generated via Equation (23) with the values L=20 
and K=19. This procedure has been developed in the MATLAB 
computer program ARMAGEN1.M. The program listing appears in 
Appendix A. 



40 



2. Double Gaussian Smoothing Function Generator 

The excitation signal chosen for generating the test 
signal response was the double Gaussian waveform, via the 
equation 

x(n) =A 1 exp(-a 1 t 2 ) -A 2 exp(-a 2 t 2 ) , ( 43 ) 



with 



where. 



3 ! = 



_ 4ln (10) 



Tl 



T x = 0 . 147 nsec 



( 44 ) 



a 2 = 



_ 4ln (10) 



Tt 



T 2 = 0 . 314 nsec 



( 45 ) 




^ 2 =^- 1 . ( 46 ) 



This waveform is a wide Gaussian pulse with a ten percent 
width of T 2 nsec subtracted from a narrow Gaussian pulse with 

a ten percent width of T t nsec. This results in a bandwidth 

of 1 to 10 GHz. Figure 12 illustrates the. double Gaussian 
waveform and Figure 13 illustrates the spectrum of the 
waveform. This procedure has been developed in the computer 
program EXCGEN.M in MATLAB. The program listing appears in 
Appendix B. 
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3. Synthetic Noise Generator 

Synthetic noise was generated to contaminate the 
signal response by adding a time series noise signal to the 
signal response. The noise was assumed to be wide sense 
stationary and white. To produce a Gaussian distribution, a 
normal distribution function was multiplied with a standard 
deviation value o , computed via the equation 



N 



o 2 = variance 



Y, [y<*> 1 



N 



SNR 

10 10 



(47) 



where, 

• y(n) is the signal response data, 

• W is the number of time points and 

• SNR is the Signal-to-Noise Ratio in dB. 

This procedure was developed in the computer program 
NOISEGEN.M in MATLAB. The program listing appears in 

Appendix C. 

4. Spectrum Estimation 

Estimation of the power spectral density (PSD) , also 
called the spectrum of the sampled signal response, is 
obtained employing the Fast Fourier Transform (FFT) . This 
spectrum may be computed via the equation 

SU)=A.|y(*) | 2 (48) 
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where , 



• N is the number of time points (power of 2) , 

• Y(k) is the FFT of the signal y(t) and 

• S(k) is the periodogram spectral estimate of y(t) . 

This procedure was developed through the computer program 
SPECTRUM. M in MATLAB. The program listing appears in 
Appendix D. 



B. SYNTHETIC SIGNAL TESTING RESULTS 

This procedure shows the weaknesses of the Cadzow-Solomon 
algorithm by using no bias compensation. Cases were examined 
for the three different SNR's of 30, 20 and 10 dB. 
Overestimation varied from 2:1 up to 5:1, e.g., for 20 true 
poles and 60 asking poles, the ratio is 3:1. The interval 
processed was composed of 200 points, starting where the 
excitation began. 

Figures 11-23 illustrates the results of this effort. 
Although the generated SNR's were 30, 20 and 10 dB, the 

resulting SNR's over the processed window were 2 dB higher. 
Observations made led to the following factors, thus improving 
results as previously obtained from the Cadzow-Solomon 
algorithm. 

• The most significant factor is to select the excitation 
starting point as the point to start the algorithm 
processing. 
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• The second significant factor is to select the system 
order or number of asking poles. In the case of 20 true 
poles, 60 asking poles gives the most accurate results for 
all SNR's employed. 

• A third significant factor is to determine the number of 
unknown zeros. For synthetic data, the position of the 
low frequency poles was obtained with better accuracy when 
this number was equal to the number of asking poles than 
with the case where the number of asking zeros was equal 
to the number of true poles. 

Figure 23 illustrates the case of K+l=20 (as compared to 
Figure 19 for K+l=60) , where (K+l) represents the number of 
asking zeros. In the case where K+l=20, more poles were 
obtained within the unit circle (as compared to the case where 
K+l=60) at the frequencies of the true poles, but with a 
different damping factor. The experimental data indicates the 
following result: 

• There is a bias of all the poles towards the unit circle 
that results in the loss of some poles. The bias is so 
large that it forces the poles to appear on the other side 
of the unit circle as noise poles. 

By examining the spectra of the synthetic signal at different 

noise levels, as illustrated in Figures 15-16, it may be 

noticed that the high frequency poles cannot be completely 

separated from the white noise spectrum when increasing the 

noise level. This is due to the fact that these high 

frequency poles carry very little energy. 

The set of equations in matrix form developed in Equation 

(29) have been developed in the computer program CADS0L1.M in 

MATLAB. The program listing appears in Appendix E. 
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Figure 11 Double-Gaussian Smoothing Function 
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Figure 13 Synthetic signal without noise 
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Figure 17 Poles in s-plane, SNR=32 dB, for synthetic signal 
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Figure 18 Poles in z-plane, SNR=32 dB, for synthetic signal 
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GigaNepers/sec +:true algorithm: *=60 unk 
Figure 19 Poles in s-plane, SNR=22 dB, for synthetic signal 
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Figure 20 Poles in z-plane, SNR=22 dB, for synthetic signal 
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GigaNepers/sec +:true algorithm: *=60 unk 

Figure 21 Poles in s-plane, SNR=12 dB, for synthetic signal 
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Figure 22 Poles in z-plane, SNR=12 dB, for synthetic signal 
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C. THIN WIRE SIGNAL TESTING 



The Cadzow-Solomon algorithm was also tested using thin 
wire measured scattering data. The results have been compared 
with the results obtained using time domain integral equation 
(TDIE) thin wire data. The poles obtained by processing the 
TDIE computed data were assumed to be correct, even though the 
program that computes the currents on the thin wire does not 
take into account the capacitance at the ends of the wire. 

1. Thin Wire Integral Equation Computed Data 

For the thin wire, the natural resonances may be 
determined by solving the integral equation that describesthe 
current flowing on the wire. The result of this type of 
simulation is the response of the wire to a specified 
excitation field. Morgan [Ref. 14] developed a time- 
domain thin wire integral equation computer routine, based on 
the formulations of Sayre and Harrington [Ref. 15]. 
The wire used for this simulation had a length of 0.1 meter 
and a radius of 1.18 mm. The back-scattering response was 
computed at three different incident aspects, ranging from 30 
to 90 degrees, with 90 degrees representing a broadside 
aspect. The excitation waveform used was the same double 
Gaussian as that used with the synthetic data, as illustrated 
in Figure 11. Figure 30 illustrates the position pole 
estimates obtained using Cadzow-Solomon algorithm and shows 
the aspect independence of the poles for three cases. The 
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five low frequency poles appeared exactly at the same position 
for all aspects. It should be noted that these results, as 
illustrated in Figure 30, appeared to be exactly the same, 
irrespective of any variations in the parameters (No. of 
points - No. of asking poles - No. of asking zeros) used in 
processing the signal. Note that for broadside illumination, 
only the odd-numbered poles appear. Because of the physical 
symmetry of both the thin wire and incident field, the even- 
numbered modes are prevented from being excited by the 
incident field. With illumination at 45 degrees aspect, a 
spectrum with adequate energy only within the bandwidth from 
1 to 8 GHz was produced, as expected. At higher frequencies, 
backscattering is suppressed because most of the energy is 
reradiated near to the specular scattering directions. 
Therefore, only the five low frequency poles are accurately 
obtained for this case. 

The TDIE program generated a response at 30 degrees 
aspect consisting of 255 points over 5 nsec, resulting in a 
sampling frequency of 50.8 GHz. In the case of 45 and 90 
degrees aspect consisting of 240 points over 5 nsec, a 
sampling frequency of 47.8 GHz resulted. These sampling 
frequencies differ from the sampling frequency of 51 GHz used 
in the measured data. 

When processing the TDIE thin wire data, points from 
160 to 200 were used to run the algorithm starting at the 
point where the excitation starts. The number of asking poles 
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ranged from 28 (for the 90 degrees aspect) to 40 (for the 30 
degrees aspect) . The number of asking zeros or feed-forward 
order of the system was either the same as the number of 
asking poles or, was calculated by determining the length of 
early-time. This early-time was computed from the formula 

t =— -(1+cosa) ( 49 ) 

e c 

where, 

• L is the length of the wire and 

• a is the aspect angle from end-on orientation. 

This value of time was then converted to the appropriate 
number of time points, as based on the sampling time interval 
of 



n b = integer 



+1 



( 50 ) 



where, 

• T is the sampling interval, and 

• n b is the feed-forward order of the system. 

This value of n b is the minimum number of asking poles that 

may be used, presenting the number of delays in the z- 
transform for the early-time of the system. 

Another consideration in the processing of these data 
was the scaling of the excitation waveform and its position 
with respect to the computed response. Although the data were 
generated using a double Gaussian waveform with a 1 Volt peak 
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amplitude, the excitation waveform had approximately the same 
peak amplitude as the response waveform. Such scaling does 
not change the frequency contents of the exciting waveform and 
gives better results in the pole extraction due to 
minimization of some of the effects of ill-conditioning in the 
data matrix. 

To position the driving waveform with respect to the 
response waveform, the time difference between the excitation 
of the first and last point of the wire may be computed. This 
time interval is represented as 

O T 

tdei ay = — -cosa, (51) 

where a is the aspect angle. The maximum absolute value of 
the response occurs when the last point of the wire is being 
excited. The information derived from the time interval, 
t del ay ' and the maximum absolute value of the response allows 

the excitation waveform to be positioned with the first point 
of the wire. 

2 . Thin Wire Measured Data 

Measurements were performed in the Transient 
Electromagnetic Scattering Lab (TESL) to test the algorithm by 
using a thin wire with the same dimensions as the one used for 
the previous simulation. A detailed explanation of the 
procedure and the techniques used for measuring the scattering 
response from the thin wire, as well as other scale model 
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targets may be found in Morgan and McDaniel [Ref. 16] 
and Bresani [Ref. 17]. 

One measurement was available for each of the aspects 
at 30, 45 and 90 degrees. The scattering response waveforms 
obtained from the measurements are illustrated in Figures 24, 
26 and 28. These waveforms were compared with the calculated 
waveforms obtained from the TDIE program. It can be seen from 
the figures that the natural modes between the calculated and 
the measured waveforms do not have exactly the same 
frequencies. 

The spectrum of the measured response waveforms were 
also obtained, giving a distribution of energy within the 
bandwidth 1-12 GHz on eight frequencies, as shown in Table 
III: 



Table III DISTRIBUTION OF ENERGY WITHIN THE SPECTRUM OF 
THIN WIRE 



Aspect 

[degrees] 


Frequencies 
[GHz] with 
most energy 


Frequencies 
[GHz] with 
less energy 


Frequencies 
[GHz] with 
no energy 


30 


3, 4. 4, 5. 8, 
7.2 


1.6,8.6,10.2 


11.6 


45 


3,4.4 


1.6, 5. 6 


7. 2, 8. 6, 10.2, 
11.6 


90 


1.6, 4. 4 


7.2 


3, 5. 6, 8. 6, 
10.2,11.6 
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Figures 25, 27 and 29 illustrate the spectrum of the measured 
thin wire scattering response for each aspect. 

When the measured thin wire response was processed, 
the points used to run the algorithm were from 180 to 200, 
starting at the initial excitation point, while the number of 
asking poles ranged from 30 to 60. The best results were 
obtained when the number of asking poles approached 40. The 
feed-forward order of the system was either the same as the 
number of asking poles, or was calculated by determining the 
early-time length. The scaling of the driving waveform and 
the positioning of the excitation with respect to the response 
was computed, as previously described. Figure 30 illustrates 
the extraction of the poles from the TDIE response, while 
Figure 31 illustrates the extraction of the poles from the 
measured waveform for the combined aspects. Figures 32-37 
illustrate the extraction of the poles from the measured 
waveform for each aspect separately. The pole results 
obtained using both the bias compensation method developed in 
this thesis and Cadzow-Solomon's bias compensation method have 
been plotted in Figures 32-37 along with the poles obtained 
without any bias compensation method, so that the results from 
the three types of methods may be compared. 
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solid : computed by TDIE dashed : measured dotted : excitation 
Figure 24 Thin wire TDIE & Measured Response, 3 0 deg aspect 




Frequency in GHz 

Figure 25 Spectrum of thin wire measured response, 30 deg 
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solid : computed by TDIE dashed : measured dotted : excitation 
Figure 2 6 Thin wire TDIE & Measured Response, 45 deg aspect 
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solid : computed by TD1E dashed : measured dotted : excitation 
Figure 28 Thin wire TDIE & Measured Response, 90 deg aspect 
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Figure 30 Poles extracted from TDIE response, all aspects 





80 




70 




60 


O 


50 


V 

J/i 




*o 




03 

tc 


40 


a 




5 


30 




20 




10 



0 . 



+ : TDIE response 

Measured response 

No bias compensation method 

o : 30 deg, x : 45 deg, * : 90 deg 



it • O 



•x-b 



-9 



-8 



-5 



-4 



-3 



-2 



-1 



GigaNepers/sec 

Figure 31 Poles from measured response, all aspects 
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Poles extracted from measured response, 45 deg 
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Poles extracted from measured response, 90 deg 
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Figure 35 Poles at z-plane, 30 deg aspect 
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Figure 36 Poles at z-plane, 45 deg aspect 
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Figure 37 Poles at z-plane, 90 deg aspect 
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IV. POLES FROM SCALE MODEL AIRCRAFT 



Scattering data for four different scale model aircraft 
targets were processed without bias compensation, using the 
Cadzow-Solomon algorithm. Poles were extracted from the 
measured scattering responses of the aircraft targets to 
double-Gaussian electromagnetic excitation incident at 0, 30, 
90, and 180 degrees. The 0 degree aspect represents a nose-on 
measurement, while the 90 degree represents a broadside 
measurement. Signatures are shown in Figures 38-45 for targets 
1 and 2 . 

The bandwidth of the TESL, 1-12 GHz, was matched to the 
scaling factor of the model aircraft targets. This scaling 
factor was 1/72 for all of the model aircraft targets used. 
Table IV contains the full scale, significant dimensions of 
these aircraft targets. The aircraft data was collected by 
Bresani [Ref. 17] at the four incident angles, as previously 
described. 



Table IV SCALED DIMENSIONS OF AIRCRAFT TARGETS 



Target number 


1 


2 


3 


4 


Overall length (cm) 


26.5 


22.5 


23.6 


26.2 


Wingspan (cm) 


19.8 


19.0 


16.4 


16.0 
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A. DATA PROCESSING 



The model aircraft scattering data was processed using a 
number of time points ranging from 200 to 240. The algorithm 
processing was started at the point of initial excitation. 

The positioning of the excitation waveform was set up 
manually with respect to the response waveform from observa- 
tions of the scattering data. Manual positioning for the 
driving waveform was required as no obvious condition existed 
for the model aircraft, as was the case for the thin-wire. 

The number of asking poles was set to 60, as the value of 
60 was found from previous experimentation to represent the 
most suitable results. The number of asking zeros was the 
same as the number of asking poles. This value was based on 
the fact that the computed early-time for the aircraft target 
is usually about 2L/c. Conversion of this early-time to the 
appropriate number of time points (equal to the number of 
asking zeros) provides a number larger than the number of 
asking poles. However, the number of asking zeros may not be 
larger than the number of asking poles, the two values are set 
to be equal. 

B. RESULTS FROM EXTRACTING THE POLES 

The model aircraft scattering data showed that the poles 
were less likely to group than the thin-wire data. The highly 
complex nature of the aircraft scatterer combined with the 
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small number of significant data points (about 150) became a 
difficult problem for the algorithm to solve. 

Figures 46-47 illustrate the spectrum of the data record 
for aircraft targets 1 and 2, for different incident angles. 

Figures 48-51 illustrate the extraction of the poles from 
the same aircraft target for the combined aspects. To obtain 
an initial indication of the possibility of target identifi- 
cation by pole extraction, nose-on measurement results have 
been compared in Figures 52-53. 

The thin wire data (Chapter III) indicated that the 
principle poles provided well defined clusters for all the 
incident angles examined, despite the use of a wide range of 
processing parameters. However, the scale model aircraft 
targets did not provide any well defined clusters under the 
limited attempts at pole estimation conducted here. The prin- 
ciple reason for this difference appears to be that the scale 
models have more complicated pole patterns than the thin wire 
targets and insufficient time was available to explore ideas 
for optimizing the performance of the estimation algorithm. 

Processing the experimental scattering data for aircraft 
models, as performed herein, is only an initial attempt to 
demonstrate resonance estimation for real-world targets 
embodying complex configurations. Processing for this 
aircraft data was conducted for only one week at the end of 
the thesis research. It is recommended that this same data be 
used for a follow-on thesis. 
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solid: response dash:excitation 

Figure 39 Measured response, model aircraft 1, 30 deg 

aspect 
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solid:response dashrexcitation 

Figure 40 Measured response, model aircraft 1, 90 deg 

aspect 




solid.response dash:excitation 

Figure 41 Measured response, model aircraft 1, 180 deg 

aspect 
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solid: response dash:excitation 

Figure 43 Measured response, model aircraft 2, 30 deg 

aspect 
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solid:response dashrexcitation 

Figure 44 Measured response, model aircraft 2, 90 deg 

aspect 




solid:response dash:excitation 

Figure 45 Measured response, model aircraft 2, 180 deg 

aspect 
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Figure 46 Spectrum of model aircraft 1 measured response, 
all aspects 
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Figure 47 Spectrum of model aircraft 2 measured response, 
all aspects 
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Figure 48 Poles extracted from model aircraft 1 measured 
response in s-plane 
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Figure 49 Poles extracted from model aircraft 1 measured 
response in z-plane 
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Figure 50 Poles extracted from model aircraft 2 measured 
response in s-plane 
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Figure 51 Poles extracted from model aircraft 2 measured 
response in z-plane 
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Figure 52 Poles extracted from four model aircraft in s- 
plane, nose-on 
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Figure 53 Poles extracted from four model aircraft in z- 
plane, nose-on 
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V. SUMMARY AMD CONCLUSIONS 



A. SUMMARY 

This thesis has attempted to demonstrate radar target 
identification by building on the earlier work of Norton [Ref. 
11], Larison [Ref. 12] and Murphy [Ref. 13]. The largest 
portion of the work consisted of testing the Cadzow-Solomon 
algorithm using synthetic and thin wire measurement data. 

Chapter I introduced the use of the Kumaresan-Tufts and 
the Cadzow-Solomon algorithms to locate target poles for a 
Non-Cooperative Target Recognition system. The resonance- 
based radar target identification problem was discussed and 
the transient electromagnetic scattering described. 

Chapter II consisted of two parts. The first part 
discussed early methods used to solve target identification 
problems. This discussion included a description of Prony's 
method and also, Singular Value Decomposition on which the 
Kumaresan-Tufts and Cadzow-Solomon algorithms were developed. 
The Kumaresan-Tufts algorithm was developed in detail, 
including Norton's bias compensation method. The performance 
of the Kumaresan-Tufts algorithm was also demonstrated, as 
illustrated in Figures 3 through 10. The second part 
described the Cadzow-Solomon algorithm. Two bias compensation 
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methods were included: that of Cadzow-Solomon and that of the 
author. 

Chapter III demonstrated the Cadzow-Solomon algorithm 
testing in two phases. The first phase of testing was 
performed with synthetic data, while the second phase was 
performed with thin wire measurement data. The thin wire data 
testing phase attempted to evaluate conditions appearing from 
a real target response. Results of the synthetic signal model 
testing are illustrated in Figures 13 through 23. The results 
of the thin wire scattering signature testing are illustrated 
in Figures 24 through 37. 

Chapter IV considered an initial attempt at extraction of 
poles from scale model aircraft measurements obtained in the 
NPS Transient Electromagnetic Scattering Lab. Measurement 
data was processed using the Cadzow-Solomon algorithm and the 
unsuccessful results are illustrated in Figures 38 through 53. 

Testing of both the Kumaresan-Tufts and the Cadzow-Solomon 
algorithms was performed using MATLAB programs. A sequence of 
programs was written to complete the demonstration of the 
target identification problem. Appendices A to E present only 
a part of this sequence of programs, including the theoretical 
models of the Cadzow-Solomon algorithm. 
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B. CONCLUSIONS 



Both the Kumaresan-Tufts and the Cadzow-Solomon algorithms 
can effectively extract poles from the scattering response of 
simple high-Q targets such as thin wires. As the late-time 
portion of the response signal is weak (with low SNR) and the 
Cadzow-Solomon algorithm has the ability to use the early-time 
portion, this thesis concentrated mainly on the use of the 
Cadzow-Solomon algorithm. 

The Cadzow-Solomon algorithm extracts poles of 
synthetically generated data, integral equation computed data, 
and thin wire measured scattering data with excellent results. 
A signal-to-noise ratio above 10 dB is required, depending on 
the damping rate of the data. The system order is 
intentionally overestimated. The excess order provides the 
flexibility to model the noise and improves the estimation of 
parameters of exponentially damped sinusoidal signals in 
noise. The Singular Value Decomposition method alleviates 
severe ill-conditioning of the data matrix. Backward 
prediction and SVD are used to separate the computed signal 
poles from spuriously computed noise poles introduced by the 
overestimated system order. 

The most critical parameters required for the successful 
thin-wire processing were the selection of the appropriate 
starting point to begin the data processing and the 
appropriate system order. The best results were obtained with 
the starting point of the data processing set at the 
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excitation starting point, and with the processing system 
order at about 3 times the true system order. 

The existence of noise in the target's response produces 
a bias in the positioning of the extracted poles. Thus, 
several bias compensation schemes may be developed. 

Given a very limited initial effort, the Cadzow-Solomon 
algorithm was unable to extract poles of scale model aircraft 
scattering data with adequate accuracy. It is conjectured 
that the data points in the natural mode information response 
are too few for the algorithm to model both the target poles 
and the noise poles. 
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APPENDIX A. ARMA COEFFICIENT AND RECURSIVE SIGNAL GENERATOR 



PROGRAM LISTING 

% ARMAGEN1.M 

% 

% Program generates a(k), b(k) coefficients of 
% 

% b(0)+b(l)*z A -l+...+b(q)*z A -q 

% H (z) = 

% 1 +a(l) *z A -l+. . .+a(p) *z A -p 

% 

% given the s-plane poles, residues, number of time points 
% and time window. 

% Programed by Gregory Lazarakos, 5 Mar 1991. 

% 

!del temp. mat 
case= 1 synt • ; 

dt=tO./ (notp-1) ; 
i=sqrt (-1) ? 

dm=exp( (sigma+i* omega) *dt) ; 
alpha=ampl . *exp (i*phase) ; 
dm=dm • ; 
alpha=alpha ' ; 

[ 'Please wait. . . ' ] 

[b,a]=residue(alpha,dm,0) ; 

A=real (a) ; 

B=real (b) ; 

% Program generates time response of an ARMA system 
% via the equation 

% N L 

% y (n) = SUM A(k) *y (n-k+1) + SUM B(k) *x(n-k+l) ,forn=l:notp 
% k=2 k=l 

% given A(k) , B(k) and input excitation record x(n) . 

% Programed by Gregory Lazarakos, 10 Mar 1991. 

N=size (A) ; 

N=N ( 2 ) ; 

L=size (B) ; 

L=L ( 2 ) ; 
for j=2:N 

A(j)— A(j) ; 

end 

% 

['Excitation is double Gaussian'] 
excgen 
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% 

[ ' Please wait. . . ' ] 



ye=zeros (l,notp) ; 
ye ( 1) =B(1) *x(l) ; 
energy= (ye (1) ) . ^2 ; 
for n=2:notp 
ye (n) =0 . 0 ; 

Ln=[n;L] ; 

Lmax=min (Ln) ; 
for k=l : Lmax 

ye(n)=ye(n)+B(k) *x(n-k+l) ; 

end 

Nn=[n;N] ; 

Nmax=min(Nn) ; 
for k=2:Nmax 

ye (n) =ye (n) +A(k) *ye (n-k+1) ; 

end 

energy =energy+ (ye (n) ) .* 2 ; 

end 

rms=sqrt (energy) ; 
ic=zeros (l,notp) ; 
for n=l:notp 

ic(n)=ye(n) ./rms; 

end 

% 

axis([0 notp -0.4 0.8]) 
plot (ic) 

title ([' Generated signal w/o noise, medium Q, 

' , int2str (nop) , ' poles')); 
xlabel('time points'); 
ylabel (' amplitude rms'); 
grid; 
pause 

hard=input( 'Do you want a hardcopy for the plot? [ 
• , • s • ) ; 

if hard== ' y ' 

plotf ile=input ( 'Enter filename : ','s'); 

eval(['meta ',plotfile]) 

end 



by 



n] : 
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APPENDIX B. DOUBLE GAUSSIAN SMOOTHING FUNCTION GENERATOR 



PROGRAM LISTING 

% EXCGEN . M 

% 

% Program generates excitation record x(n) for a response 
% from TESL data and synthetic data, 

% via the equation of a double Gaussian waveform 
% 

% x(n) = Al*exp(-al*t A 2) - A2*exp(-a2*t A 2) 

% 

% given the time window. 

% Program uses input values of the 10% height of the low 

% and high frequency ends of the double Gaussian frequency 

% response to determine the pulse widths in the time domain. 
% Programed by Gregory Lazarakos, 8 Apr 1991. 

% Last Revision August 6, 1991. 
if case=='meas' 

if filename (1)=='F' 

LENGTH=0 .25; 

elseif filename(siz-l:siz)=='sp' 

LENGTH=eval (filename (3 : 4 ) ) ; 

LENGTH=LENGTH/ 100 ; 

else 

LENGTH=0 . 1 ; 

end 

end 

npts=notp ; 
tmin=0 . 0 ; 
tmax=t0 ; 

DT=dt ; 

T1=0 .147 ; 

T2=0 .314; 

thr=input( ' Enter threshold value for time in nsec 
[1.25] : ' ) ; 

if isempty(thr) 
thr=1.25; 

end 

al= (4 . 0 . *log(10) )./(Tl. A 2); 
a2= (4 . 0 . *log(10) ) ./ (T2 . A 2 ) ; 

Al=sqrt(al) ./(sqrt(al)-sqrt(a2) ) ; 

A2=A1-1 . 0 ; 
if case=='synt' 

point=input ( 'Which to be the point, the excitation 
starts? ' ) ; 

x=zeros(l,npts) ; 
for n=l:npts 
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t= (n-ll-point) . *DT; 
if t<=thr 

el=Al. *exp (-al . *t . A 2) ; 
e2=A2 . *exp (-a2 . *t. A 2) ; 
x(n) =el-e2 ; 
else 

x(n) =0.0; 

end 

end 

elseif case=='meas' 

asprad=aspect* (pi/180) ; 

top=input ( 'Do you want to position the excitation, auto or 
manual ? (a,m) : ','s'); 

if top=='a' 

[maxic, indexl]=max(ic) ; 

[minic, index2]=min(ic) ; 
x=zeros(l,npts) ; 

tdel=cos (asprad) * (2*LENGTH/3e8) ; 
ndel=round (tdel/ (dt*le-9) ) ; 
if abs (maxic) >=abs (minic) 
point=indexl-ll-ndel ; 
for n=l:npts 

t= (n-indexl+ndel) .*DT; 
if abs (t) <=10*DT 

el=Al. *exp(-al . *t. A 2) ; 
e2=A2.*exp(-a2.*t. A 2) ; 
x(n)=el-e2; 
else 

x(n)=0.0; 

end 

end 

x=x* 0.6; 
else 

point=index2-ll-ndel ; 
for n=l:npts 

t= (n-index2+ndel) .*DT; 
if abs (t) <=10*DT 

el=Al.*exp(-al.*t. A 2) ; 
e2=A2 . *exp (-a2 . *t. A 2) ; 
x(n)=el-e2; 
else 

x(n)=0.0; 

end 

end 

x=x*0.6; 

end 

end 

if top=='m' 

point=input ( 'Which to be the point, the excitation 
starts ? ' ) ; 

x=zeros (l,npts) ; 
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for n=l:npts 

t= (n-ll-point) . *DT; 
if abs (t) <=10*DT 

el=Al. *exp(-al. *t . A 2) ; 
e2=A2 . *exp (-a2 . *t . A 2) ; 
x(n)=el-e2; 
else 

x(n)=0.0; 

end 

end 

x=x* . 6 ; 

end 

ltst= ( l+cos (asprad) ) * (LENGTH/3e8 ) ; 
nltst=round(ltst/ (dt*le-9) ) ; 
ltimp=point+22+nltst ; 
else 
exit 
end 

if case=='synt' 

axis([0 notp -0.4 1]) 
plot (x) 

title (' Double-Gaussian Smoothing Function'); 
xlabel('time points'); 
ylabel ( 'amplitude rms'); 
elseif case=='meas' 

plot (1 :npts, x, ' — ' , l:notp, ic, It imp, ic (ltimp) , ' * ' ) 
if filename (1)=='F' 

title ( [ ' Data response signal from aircraft' , filename (1 : 3 ) , ' , 
aspect= ', num2str (aspect) , ' deg ']); 

elseif filename(siz-l:siz)=='sp' 
title ( [ ' Data response signal from sphere ' , filename (3:4) , 'em, 
aspect*' ,num2str (aspect) , ' deg ')); 
else 

title(['Data response signal from thin wire, 
aspect*' ,num2str (aspect) , ' deg ']); 
end 

xlabel ( ' solid=response dash=excitation ' ) ; 

ylabel ( 'amplitude '); 

text ( ltimp, 0. 1, ' late time >'}» 

if filename(l:2)=='et ' 

text (100, 0 . 4 , 'TDIE response'); 
end 
grid ; 

else 

exit 

end 

pause 

hard=input ( ' Do you want a hardcopy for the plot? 
[n] : ' , 's' ) ; 

if hard** ' y ' 
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plotf ile=input ( 1 Enter filename : 
eval(['meta plotf ile]) 



• , ' s ' ) ; 



end 

if case=='meas ' 

save temp ic energy x point dt to notp filename aspect siz 
LENGTH 
end 
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APPENDIX C. SYNTHETIC NOISE GENERATOR 



PROGRAM LISTING 

% NOISEGEN.M 

% 

% Synthetic Noise Generation 
% Programmed by Greg Lazarakos 20/2/91 
% 

SN=input ( ' How many db ? (7,10,20,30) [20]:'); 

if isempty(SN) 

SN=20; 

end 

sn=SN./10; 
sn=10. A sn; 
en=0 ; 

for k=l:notp 

en=en+ ( ic (k) ) . A 2 ; 

end 

ava=en./ (notp. *sn) ; 
ava=sgrt (ava) ; 
rand ( 'normal ' ) 
w=rand ( ic) *ava ; 
ic=ic+w; 

% 

save temp ic x point filename w ye energy dt dm notp omega 
sigma to 
% 

axis([0 notp -0.4 0.8]) 
plot (ic) 

title ( [ 'Generated signal with noise S/N=' , int2str(SN) , ' db' ] ) ; 
xlabel('time points'); 
ylabel ( ' amplitude rms ' ) ; 
grid ; 
pause 

hard=input ( ' Do you want a hardcopy for the plot? 
[n]: ', ’s') ? 

if hard== ' y ' 

plotf ile=input (' Enter filename : ','s'); 

eval ( [ ' meta ' , plot f ile ] ) 

end 
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*> <N> <H> <H> 



APPENDIX D. SPECTRUM ESTIMATION 



PROGRAM LISTING 

SPECTRUM. M 

Programed by Gregory Lazarakos, 16 Apr 1991. 

Last revision, 6 Aug 1991. 
df=l ./ ( (notp-1) . *dt) ; 

IC=fft (ic,notp) ; 
f=(0: notp-1) *df ; 
fmax=l./dt ; 
axis([0 fmax -3 3]) 
plot (f ,IC) 

title (' Fourier transform of the response signal'); 
xlabel ( ' Frequency in GHz ' ) ; 
y label ( 'Amplitude ' ) ; 
grid ; 
pause 

SIC= (abs (IC) ) . A 2 ; 

S I C=S I C/ no tp ; 
axis ( [0 20 0 max (SIC) ] ) 
plot(f ,SIC) 
if isempty(SN) 

title ( 'Spectrum of the response signal '); 
else 

title ([' Spectrum of the response signal 
S/N= ' , int2str (SN) , ' db')); 
end 

xlabel (' Frequency in GHz'); 
y label ( ' Amplitude ' ) ; 
grid ; 
pause 

hard=input ( ' Do you want a hardcopy for the plot? 
[n] : ' , ' s ' ) ; 

if hard== ' y ' 

plotf ile=input (' Enter filename : ','s'); 

eval ( [ ' meta ' , plotf ile ] ) 

end 



with 
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APPENDIX E. THE CADZOW-SOLOMON POLE EXTRACTION ALGORITHM 



PROGRAM LISTING 

% CADSOL1.M 

% 

% Cadzow-Solomon' s Algorithm for Extracting the Poles 
% and Residues 
% Version 1.0 

% Programmed by Greg Lazarakos 3/16/91 
% 

% Set the first sample (kapa) 
if point>0 

kapa=point ; 
else 

kapa=l ; 

end 

% Set the number of samples (CN) 

CN=input ( 'How many samples to process ? (>120) [200]:*); 

if isempty(CN) 

CN=200; 

end 

% Compute the SNR for the processed time window 
en=0 ; 
nen=0 ; 

for k=kapa:CN 

en=en+(ic(k) ) . A 2; 
nen=nen+ (w(k) ) . A 2 ; 

end 

SNR=en./nen; 

SNR=10 . *logl0 (SNR) ; 

% Set the number of poles (CL) > number of real poles (CM) 
% CM <= CL and 2*CL <= CN-CL 

CL=input ( 'How many unknown poles ? (>20) [60]: '); 

if is empty (CL) 

CL=60; 

end 

% Set the number of expected zeros (CK) 

CK=input ( 'How many expected zeros ? [No. poles]: '); 

if isempty(CK) 

CK=CL 

end 

['Please wait...'] 

% Forming the matrix YI [ (CN-CL) *CL] from the data ic(n) 
YI=zeros (CN-CL, CL) ; 
for j=kapa+l:CN-CL+kapa 

YI ( j-kapa, : )=ic( j : j+CL-1) ; 

end 

% Forming the matrix h[ (CN-CL)*!] from the data ic(n) 
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h=[ ] ; 

h=ic (kapa : CN-CL+kapa-1) ; 
h=h • ; 

% Forming the matrix XI [ (CN-CL) *CK] from the data x(n) 
XI=zeros (CN-CL, CK) ; 
for j=l: CN-CL 

XI ( j / : ) = x ( j : j+CK-1) ; 

end 

% Unify matrices YI and XI as I=[Yl|XI] 

I=[YI XI] ; 

% Set the tolerance 
tol=0. 000001; 

% Find the vector of backward prediction coefficients 
beta=pinv(I, tol) *[h] ; 

% Set the coefficients of the prediction-error filter 
polynomial 

ca=zeros(l,CL+l) ; 
ca (1)=1; 
for j=l : CL 

ca ( j+l)=-beta ( j ) ; 

end 

% Set the coefficients of the polynomial B(z) 
cb=zeros(l,CK) ; 
for j=l:CK 

cb ( j ) =beta ( j +CL) ; 

end 

% Find the residues and the poles 
[resid / d,gk]=residue(cb,ca) ; 

% Find the signal poles 
s=log (d) /dt; 
k=0 ; 

1 = 0 ; 

for j=l:CL 

if real (s ( j ) ) >0 , 
k=k+l; 

polesl(k)=-s(j) ; 
residl(k)=resid(j) ; 
else 

1 = 1 + 1 ; 

nopoll(l)=-s(j); 
noresl(l)=resid(j) ; 

end 

end 

dml=exp (polesl*dt) ; 
dm2=exp (nopoll*dt) ; 
alphal=exp(residl*dt) ; 

CM=k; 

[ 'The signal has the following ' ,num2str(CM) , ' real poles. ' ] 
polesl 
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